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Abstract. We describe a general framework for avoiding spurious eigenvalues — unphysical 
unstable eigenvalues that often occur in hydrodynamic stability problems. In two example prob- 
lems, we show that when system stability is analyzed numerically using descriptor notation, spu- 
rious eigenvalues are eliminated. Descriptor notation is a generalized eigenvalue formulation for 
differential-algebraic equations that explicitly retains algebraic constraints. We propose that spu- 
rious eigenvalues are likely to occur when algebraic constraints are used to analytically reduce the 
number of independent variables in a differential-algebraic system of equations before the system is 
approximated numerically. In contrast, the simple and easily generalizable descriptor framework si- 
multaneously solves the differential equations and algebraic constraints and is well-suited to stability 
analysis in these systems. 
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1. Introduction. Spurious eigenvalues are unphysical, numerically-computed 
eigenvalues with large positive real parts that often occur in hydrodynamic stability 
problems. We propose that these unphysical eigenvalues occur when the incompress- 
ible Navier Stokes equations are analytically reduced - i.e., the algebraic constraints 
are used to reduce the number of independent variables before the system is approx- 
imated using spectral methods. 

An alternative approach to analyzing differential-algebraic equations is the de- 
scriptor framework, posed as a generalized eigenvalue problem, which explicitly re- 
tains the algebraic constraints during the numerical computation of eigenvalues. We 
reformulate two common hydrodynamic stability problems using descriptor notation 
and show that this method of computation avoids the spurious eigenvalues gener- 
ated by other methods. The descriptor formulation is a simple, robust framework 
for eliminating spurious eigenvalues that occur in hydrodynamic stability analysis. 
Additionally, this formulation reduces the order of the numerically approximated dif- 
ferential operators and accommodates complex boundary conditions(BCs), such as a 
fluid interacting with a flexible wall. 

Resolving the spectrum of hydrodynamic operators is critical for time integration, 
linear stability [8] and transient growth analyses [130I13II3II1ISII3IIBII6] of 
fluid flows. Spurious eigenvalues often arise in these spectral numerical computations 
describing incompressible fluids and must be identified or eliminated. Although the 
term "spurious" is applied to many numerical anomalies, it is used here to refer ex- 
clusively to large, positive eigenvalues that arise due to application of extra boundary 
conditions. These are distinct, for example, from spurious pressure modes that oc- 
cur when the momentum equation is collocated at interior Chebyshev-Gauss-Lobatto 
nodes [28] , 

Researchers have developed special methods to avoid or filter these modes and 
uncover the true spectrum of the model problem. Perhaps the first description of 
these unphysical values is given by Gottlieb and Orszag [13]. Many other researchers 
have encountered similar modes [U [32J, [39] and developed methods for avoiding [TBI 
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[TTllL 25, 38 or filtering [57112 them. These methods for avoiding spurious modes 
are specific to very special clamped BCs where homogeneous Dirichlet and Neumann 
conditions hold at the boundaries. We emphasize that the descriptor method can be 
implemented using widely available numerical libraries and is easily generalized to 
non-clamped BCs. 

This paper is divided as follows: Section [5] discusses the general descriptor frame- 
work and Sections [3] and 2] develop this method for specific examples, the Orr- 
Sommerfeld operator and Gottlieb and Orszag's one-dimensional potential flow model [13] . 
Section [5] discusses how infinite eigenvalues arise in formulations of the incompress- 
ible fluid equations, and explains how the descriptor formulation explicitly accounts 
for these eigenvalues. We also discuss why methods that include analytical transfor- 
mations might generate spurious eigenvalues, as well as benefits and drawbacks of 
descriptor formulations. 

2. The descriptor framework. Descriptor notation can be used to describe 
any dynamical system where a set of differential-algebraic equations is reduced to a set 
of differential equations only. In this paper we focus on two models for incompressible 
fluid flow: the 2D Navier-Stokes equations which generate the Orr-Sommerfeld (OS) 
operator, and Gottlieb and Orszag's ID potential flow model. 

Descriptor notation was developed in the control theory community J24[[3TJ[22] to 
describe and analyze systems of differential-algebraic equations. In descriptor form, 
the differential time operator is preceded by a square, possibly singular matrix: 

E^- t 4> = A ( f>. (2.1) 

In descriptor systems, stability is determined by the generalized eigenvalues of the 
(A,E) system, which is the ratio of the pair (a, (3) where (3Au = aEu for some 
non-zero vector u. While traditional eigenvalues are never infinite, descriptor system 
can have (many) infinite eigenvalues. If E contains a zero row corresponding to an 
algebraic constraint, there will be an infinite eigenvalue corresponding to the infinitely 
fast dynamics of that constraint. 

Let us assume that we have a system of differential-algebraic equations for n 
fields. Let there be m algebraic constraints, and k = n — m equations that contain a 
differential time operator. Physical systems are often modeled by differential-algebraic 
equations of this form because algebraic constraints often arise as approximations to 
differential equations for quickly equilibrating variables. 

Let the vector v contain the k fields that are acted upon by the differential time 
operator, and p contain the remaining m fields. Each of the n fields can be discretized 
using N points for each field, which results in a finite-dimensional system of differential 
algebraic equations that can be written as follows: 



I 


E4> = A<j>. (2.3) 

The discretized operators A and E are square matrices with (n x iV) 2 entries. 
If the original equations are partial differential equations, the operator A contains 
spatial derivatives that require BCs. Any BC can be incorporated into the numerical 
solution using boundary bordering in a straightforward manner [3J. For example, if 
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A contains spatial derivatives up to order G for each of the k fields in v, then 67 BCs 
are required for each field and are enforced by using (67 x k) rows of A as additional 
algebraic constraints. We then solve for the generalized eigenvalues of the (A, E) pair. 
We assume the pair is regular, which means that det(sE — A) is not identically zero for 
all s. There are many numerical routines that solve for the generalized eigenvalues of 
regular matrix pairs when one of the matrices is singular (see references in [36] and [1]). 
A widely available routine is the MATLAB (LAPACK) QZ' algorithm [26]. Note that 
E contains (mN) zero rows corresponding to the original algebraic constraints and 
(kG) zero rows corresponding to the boundary bordering constraints, resulting in an 
(A, E) pair with (mN + kG) infinite generalized eigenvalues and (kN — kG) finite 
eigenvalues. Although the algebraic constraints are approximated using boundary 
bordering, descriptor notation ensures that the generalized eigenvalues corresponding 
to these constraints are not approximated - they are formally infinite. 

Our method is in contrast to other methods where the algebraic constraints are 
removed analytically. In these methods, the system is analytically converted to a 
system of equations that contains only fields in v, and consequently, a subset of the 
algebraic constraints on p are converted to additional BCs on v. Discretization and 
enforcement of the new BCs may lead to an approximate enforcement of the original 
algebraic constraints, resulting in spurious eigenvalues. 

One example of this over-specification of the BCs occurs in the derivation of 
the OS operator. Assume that A±% contains a first order spatial derivative and 
A22 is zero. We eliminate p, using a method described for the OS operator in Ap- 
pendix [A] (Eqs. IA.UA.7|) . and the resulting system of equations contains a spatial 
derivative of order 2 + 67 acting on the fields in v. Because the derivative operator is 
two orders higher than before, the system requires two new BCs. These extra BCs are 
determined by numerically approximating the algebraic constraints at the boundary. 
However, the algebraic constraints were used to eliminate p and were already eval- 
uated at the boundary in the analytical computation. In several spectral methods, 
the algebraic constraints at the boundary are enforced twice, with slightly different 
numerical approximations in each case. 

We now illustrate how the descriptor method avoids this over-specification in two 
example problems. 

3. Example: Orr-Sommerfeld Operator. The Navier-Stokes equations for a 
compressible, viscous fluid can be written [23] : 



where v represents the three components of the velocity, p is the fluid pressure, p is 
the the mass density, and v and £ are the first and second kinematic viscosities, re- 
spectively. An additional thermodynamic constraint such as constant entropy relates 
the variables v, p and p and closes the system of equations. 

As the system approaches the incompressible limit, the partial differential equa- 
tion expressing conservation of mass (Eq. 13. 2\ becomes increasingly stiff, and in the 
limit a differential equation in time is replaced by an algebraic constraint, 




(3.1) 



(3.2) 



y-v = o. 



(3.3) 
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As the equation for the mass density 13.21 becomes more stiff the eigenvalues of the 
operator will have larger real parts, and in the limit the eigenvalues will be infinite. In 
this case the pressure can be thought of as a Lagrange multiplier that instantaneously 
satisfies the divergence constraint and no thermodynamic equation is required to 
relate the pressure and density. Therefore Eqs. (|3.1[ 13. 3[) alone are the Navier-Stokes 
equations for an incompressible fluid. 

Numerical solutions to the linearized, incompressible equations of motion have 
traditionally been determined by combining the Laplacian and divergence of the NS 
equations. This allows the pressure to be rewritten in terms of the fluid velocities, 
thereby reducing the set of four differential algebraic equations to two differential 
equations. For later reference, these equations are derived using operator multiplica- 
tion in Appendix [A] In planar channel flow this results in the OS equation for the 
wall-normal velocity and the Squire equation for the wall- normal vorticity |211 115] . 
Descriptor notation does not require this reduction, and does not generate the fourth 
order OS operator. 

In most two-dimensional flow models the system is assumed to be translation 
invariant in the spanwise(z) direction, so the coupling between the vorticity and ve- 
locity is zero and eigenvalues of the OS operator determine the system stability. In 
order to compare our results with previous studies, we focus here on the stability of 
the OS operator alone. 

The OS operator contains the term A 2 , which is a fourth-order spatial derivative 
in the wall-normal direction and requires four BCs. Two BCs are simply the no-slip 
conditions from the original equations, v y (±l) = 0. The remaining two BCs arise 
from the divergence constraint: 

J=0 + 0, (3.4) 

where the last equality holds because v x and v z are constant (zero) in the streamwise 
(x)- and spanwise(z)- directions at the boundary. Therefore, homogeneous Neumann 
BC on the wall-normal velocity is a direct consequence of the incompressible limit. 
This has important implications for numerical approximations, as well will discuss in 
Section [5l 

Spurious eigenvalues arise frequently in the analysis of this OS operator I12[ 
[39] . While there are several methods for avoiding or filtering these eigenvalues, such 
as applying only Dirichlet conditions to the second order operator [37J [T7J , they 
are tailored to clamped boundary conditions and are not easily generalizable. Our 
approach to solving the original system of Eqs. (|3.11 13.3|) is more general. We avoid 
combining the two equations into a single system and instead write the system using 
descriptor notation. A similar analysis using descriptor notation was applied previ- 
ously to incompressible Stokes flow in the context of systems control 33, 34J, which 
we generalize to eliminate spurious eigenvalues. The NS equations (|3.1[ 13. 3[) can be 
written: 

I 0' 


E<j) = ~A<t>; (3.6) 
v y (±l) = » x (±l) = v z (±l) = 0. (3.7) 

The operators A, Q and D are defined in Appendix [A] 
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Numerical Method 


Ai 


A 2 


A 3 


Descriptor Chebyshcv Tau 


+.0037 


-0.0348 


-0.0350 


OS Chcbyshev Tau 


97.557 


85.735 


0.037 



Table 3.1 

This table compares the first three non- infinite eigenvalues for the 2D NS equations with a = 1, 
R = 10 000 for a descriptor vs. an Orr Sommerfeld(OS) spectral collocation scheme with N = 34- 
The OS scheme generates two large, positive eigenvalues, while the descriptor scheme does not. 



In hydrodynamic channel flow, a Fourier transform can be taken in the translation 
invariant streamwise and spanwise directions. The velocities have Dirichlct BCs at 
the channel boundaries, y — ±1, and we discretize the system of equations in the 
wall normal direction using Chebyshev collocation or other spectral methods. No-slip 
BCs are enforced for the velocities, but care must be taken to ensure that no BCs are 
applied to the pressure. 

For simplicity, boundary bordering is used to enforce Dirichlet BCs on each com- 
ponent of the velocity. We approximate each component of the velocity and the 
pressure by a vector of TV points. The BCs require that the first and last entry of 
each velocity vector is zero: 

«i(yo) =«<(+!)= 0; (3.8) 
Vi{V(N-i)) =Vi(-l) = 0; (3.9) 
i = x,y,z. 

This is equivalent to deleting the first and last columns of Cheb^ 1 ' and Cheb*- 2 ' that 
occur in the operators D and A, respectively, in Eq. 13.51 Additionally, time derivatives 
of the velocity evaluated at the boundary are zero Vi(yo) — i)i(yjsr-i) = 0, which is 
equivalent to deleting the corresponding rows of A and Q in Eq. 13.51 As a result, 
each component of the velocity is represented by an N — 2 column vector. Let M = 
3 x (N — 2). Then the operator A is approximated by an M-by-M matrix, Q is 
a matrix that is M rows by N columns, D is N rows by M columns. Note that 
there are no BCs imposed on p, which is still a vector of size N. We then solve for 
the generalized eigenvalues of the (A, E) pair using the MATLAB (LAPACK) QZ' 
algorithm [26 . In the parameter range studied, the pair is regular, so that det(s_E — A) 
is not identically zero for all s. Because E contains N rows that are singular, there 
are N infinite generalized eigenvalues and M finite eigenvalues. 

Table [3] 1 lists the first three non-infinite eigenvalues calculated from the descrip- 
tor system given by Eqs. I3.5H3.7I with 34 collocation points, a = 1 and R = 10 000. 
These three eigenvalues match those calculated using other methods which have been 
identified as real, physical eigenmodes - there are no spurious eigenvalues. For com- 
parison, Table [31 1 also lists the first three non- infinite eigenvalues generated using a 
traditional Orr-Sommcrfcld Chebyshev Tau method. The traditional method gener- 
ates a pair of large, positive, spurious eigenvalues. 

A descriptor formulation can be used for any system of differential algebraic equa- 
tions. As another example we now formulate Gottlieb and Orszag's one-dimensional 
potential flow model using descriptor notation and show that this method avoids 
spurious eigenvalues. 

4. Example: Incompressible limit of a ID fluid model. Spurious eigen- 
values have been studied in greatest depth using the model problem of Gottlieb and 
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Orszag |13j . This is a model for a two-component, one-dimensional fluid How at low 
Reynolds number, and is described by the following equations: 

Here £ is the vorticity and ijj is the stream function defined by (v Xl v y ) = (—dtp/dy, dip/dx). 
The divergence constraint is automatically satisfied because V • (v x ,v y ) = equates 
the mixed partial derivatives of ip, which is always true for analytic ip. For a fluid 
between stationary rigid walls, no-slip conditions on the velocity at the boundary 
correspond to the following constraints on the stream function: 

tp(x = ±1, t) = ip x (x = ±1, t) = 0. (4.3) 

The usual method used to find the eigenvalues of Eqs. |4~T1 combines the two equations 
into a single partial differential equation, 

Ipxxt = Wpxxxx, (4-4) 

then inverts the second order differential operator and represents the operators in 
terms of spectral differentiation matrices, 

ipt = ^{d X x)~ 1 d xxxx ip. (4.5) 

The second order operator is rendered invertible through the application of two BCs. 
However, there are four BCs on tp, so there is ambiguity as to which two should 
be applied. One choice uses basis recombination so that each of the basis functions 
individually satisfies all four BCs. Dawkins, Dunbar and Douglass [7] have shown 
that a numerical scheme with Chebyshev basis functions generates large spurious 
eigenvalues, while a scheme with Legendre basis functions generate formally infinite 
eigenvalues. They also show that Neumann BCs exactly match the form of the Legen- 
dre polynomials, and conclude that spurious eigenvalues are approximations to these 
infinite eigenvalues that occur when a Chebyshev basis is used instead of a Legendre 
basis. 

Because the BCs are algebraic constraints that correspond to infinitely fast pres- 
sure modes, the Legendre polynomials are an exact basis for approximating the al- 
gebraic constraints, and therefore exactly recover the infinite eigenvalues that corre- 
spond to these constraints. The Legendre basis recombination method is similar to 
the descriptor framework in that the infinite eigenvalues are exactly computed in each 
case. 

An alternate method for incorporating BCs is the traditional Chebyshev- Tau 
method, where boundary bordering is used to replace four terms in the Chebyshev 
expansion with four algebraic constraints. These algebraic constraints are then used 
to reduce the total number of equations by four, and the eigenvalues of the result- 
ing system of equations are computed [llj . This reduced set of differential equations 
has been shown to be equivalent to the Chebyshev basis recombination method de- 
scribed above [7], and generates spurious eigenvalues, as shown in Table 01 1. Various 
other approaches to solving the spurious eigenvalue problem involve imposing only the 
Dirichlet BCs on the second order differential operator |T7] . One such approach 
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is Weideman and Reddy's [37] spectral collocation scheme with modified clamped 
conditions, which is also listed in Table 2)1. 

In contrast to these approaches, the system Eqs. (|4.1l 14. 2p can be written in 
descriptor form: 



I 0" 




"C" 




vd xx 




"C 







J' 




—I 8 

1 u xx 







E 



"C" 



= A 



"c" 



(4.6) 
(4.7) 



A second order operator acts on ip and there are four BCs for ip, while a second order 
operator acts on £ and there are no BCs on £. This is a similar situation to the OS 
problem because the four conditions on ip come from a combination of the no-slip 
condition and the divergence constraint. 

Therefore we expect these four modes to have infinite eigenvalues and impose all 
the BCs on ip as algebraic constraints. The resulting approximation to the bottom 
row of the matrix Eq. 14.61 (0 = — £ + d xx ip) is: 

= ip ; (4.8) 
= -SijCj + Chebg } ipj,j, i S (1, N - 2); (4.9) 
= V(tf-x), (4.10) 

and the approximation to the top row of Eq. 14.61 dtC = dxxCi 1S: 

= C3hebg ) Vj, je(l,N-2); (4.11) 
8td= Cheb^O, ie(l,JV-2),iG(0,JV-l); (4.12) 
= Cheb$_ 2)j ^ j€(l,N-2). (4.13) 

With descriptor notation there are (N — 2) infinite eigenvalues corresponding to the 
(N — 2) interior points of the algebraic constraint Eq. 14.91 The four new BCs cor- 
respond to four formally infinite eigenvalues in the numerical spectrum. This final 
descriptor system contains an E matrix with (N — 2) + 4 rows of zeros, and therefore 
the system has N + 2 generalized infinite eigenvalues. If a Chebyshev-Tau method is 
used to discretize the differential operators, descriptor notation is equivalent to the 
method suggested by Gottlieb and Orszag to avoid spurious eigenvalues when they 
first posed this simple model [13j . except that we use matrix methods for singular 
matrices instead of a shooting algorithm to determine eigenvalues. Again, the 'QZ 1 
routine in MATLAB is used to compute the spectrum for this system of equations. 

Table SJl compares the first three non-infinite numerical eigenvalues to the analyt- 
ically computed eigenvalues for several different discretization schemes, and confirms 
that this method generates no spurious eigenvalues. 

Figure [4~T1 shows the eigenvalues of the Gottlieb-Orszag model for 20 discretiza- 
tion points (Fig. I4.lf a)) and 40 discretization points (Fig. 14.1T b)). and confirms that 
the descriptor method eliminates spurious eigenvalues. Open circles, computed using 
a descriptor method, are all in the left-hand plane. In contrast, the points com- 
puted using the general Chebyshev tau method (closed circles) include large spurious 
eigenvalues that increase with increasing TV. 
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Numerical Method 


Ai 


A 2 


A 3 


Exact 


-9.8696 


-20.1907 


-39.4784 


Descriptor Chebyshev Tau (collocation) 


-9.8690 


-20.1883 


-39.4694 


Descriptor Chebyshev Tau (basis fn) 


-9.8696 


-20.1907 


-39.4784 


Traditional Chebyshev Tau (basis fn) 


56,119 


48,515 


-9.8696 


Modified clamped Chebyshev Galcrkin (collocation) [37 


-9.8696 


-20.1907 


-39.4784 



Table 4.1 



Comparison of the first three non-infinite eigenvalues for the ID viscous fluid model for several 
discretization schemes with N = 20. Both descriptor schemes (collocation and basis function) 
eliminate spurious eigenvalues, while the traditional Chebyshev Tau method includes large spurious 
eigenvalues in the right-hand plane. Weideman and Reddy \31\ have developed a Galerkin collocation 
scheme that avoids spurious eigenvalues by using modified clamped boundary conditions; this method 
is specific to homogeneous BCs. See AnnendixlB] for definitions of Tau vs. Galerkin schemes and 
collocation vs. basis function schemes. 



N = 20 




-2 2 4 

Re(scaled eigenvalues) 



Fig. 4.1. Eigenvalues, scaled by x — > sign(x) + log 10 (l + |x|), of the Gottlieb- Orszag model 
for N = 20(a) and N = 40 (b) discretization points. Open circles are computed using a descriptor 
Chebyshev Tau basis function scheme, while closed circles are computed using a traditional Cheby- 
shev Tau basis function scheme. Two large, positive spurious eigenvalues occur in the traditional 
scheme and increase with increasing N . They are eliminated in the descriptor scheme. 



5. Conclusions. The descriptor framework is a generalized eigenvalue method 
for hydrodynamic stability problems that eliminates spurious eigenvalues, as shown 
in Tables [3] and [4j It ensures that infinitely fast modes will retain formally infinite 
eigenvalues, even when those eigenvalues are computed numerically. In incompressible 
fluid flows, there are N discretized pressure variables which correspond to TV gener- 
alized eigenvalues which are explicitly infinite. This method enforces only the no-slip 
BCs and applies them in an intuitive and unambiguous way. 

While other methods reduce the number of fields (from four to two in the Orr- 
Sommerfeld formulation), they do so at the expense of creating higher order deriva- 
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tive operators, which decreases resolution for a fixed number of grid points or basis 
functions. In addition, the descriptor formulation does not require inversion of a dif- 
ferential operator, is adaptable to different discretization schemes, and can be simply 
extended to problems where BCs are non-trivial. 

The descriptor framework and the traditional Orr-Sommerfeld method with mod- 
ified clamped boundary conditions are complementary schemes. The latter is likely to 
be useful for problems with homogeneous boundary conditions where the number of 
discretization points is very large. The former is more flexible: it is useful when higher 
order operators are undesirable, inhomogeneous boundary conditions are required, or 
direct calculations for the pressure (or any Lagrange multiplier in a general DAE) are 
desired. Also, newer versions of the QZ algorithm are quite efficient [SHU], making 
computations with large matrices more feasible. 

One area of research where descriptor notation is extremely promising is the study 
of a fluid interacting with a compliant boundary. In this system, no-slip conditions at 
the wall require that the fluid velocity match the wall velocity there. These compli- 
cated Dirichlet BCs can be applied directly to the second order differential operators 
in A. Furthermore, the fluid pressure at the boundary remains as an independent 
variable in the eigenvalue computation, which is advantageous because the pressure 
at the boundary influences wall motion. This topic is currently under investigation. 

Utilizing descriptor notation is similar in spirit to several other methods for avoid- 
ing spurious eigenvalues which require the algebraic constraints to be discretized sepa- 
rately from the differential equations [T3l [TTl [3] . Gottlieb's method utilizes a shooting 
algorithm for determining eigenvalues - his algorithm was developed before efficient 
matrix methods were available for solving generalized eigenvalue problems with singu- 
lar matrices. Gardner and Boyd also describe methods where the algebraic constraint 
is discretized. The descriptor framework generalizes these ideas and presents a simple, 
systematic method for avoiding unphysical spurious modes by using new and efficient 
'QZ' algorithms. Recently, Fornberg [10] has developed a fictitious point method 
for avoiding spurious eigenvalues that is applicable to problems with inhomogeneous 
boundary conditions. 

Although we have focused here on hydrodynamic stability problems, descriptor 
notation might be advantageous to any researchers who study stability of differential- 
algebraic equations using spectral methods. The simplicity and generalizability of 
the descriptor framework suggest that it is well-suited to stability analysis in many 
differential-algebraic systems, including but not limited to incompressible fluids. 
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Appendix A. Derivation of the Orr-Sommerfeld operator using matrix 
multiplication. For a channel flow between rigid walls, the nondimensionalized, 
linearized Navier-Stokes equations can be written schematically as: 



where Q is the column operator — {d x , d y , d z }', and V is the row operator {d x , d y ,d z }. 
A no-slip condition at the boundaries requires v = 0, while there are no explicit BCs 



v = Av + Qp; 
Vv = 0, 



(A.l) 
(A.2) 
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on p. The operator A can be written as 



follows: 



.4 



u 

| + Ud x 









(A.3) 



where R is the Reynolds number and U is the mean flow. The mean flow is in the 
rr-direction, while the channel walls ensure the flow is non-periodic in the y-direction. 
To derive the OS equation, we first to rewrite the pressure in terms of the velocities. 
Because the operator T> commutes with time derivatives, left multiplication by T> on 
Eq. I A. 1 1 results in a left-hand side which is identically zero. Therefore p and v have 
the following relationship: 



V {Av) = -VQp. 



(A.4) 



A three-by-three matrix is formed from the product of the scalar VQ = —A = 
— (d x + dy + dl) and the identity matrix (7). Left multiplication by AI of Eq. IA.ll 
results in the following equation: 



AIv = AI ( Av) + AIQp. 



(A.5) 



Because AI commutes with Q, the two may be interchanged in the last term of 
Eq. IA.51 The pressure p can then be removed from the equation using Eq. IA.4I 
resulting in: 



AIv = AI (Av) + (QV) (Av) , 



(A.6) 



where QT> is a three-by-three matrix. Rewriting the velocity fields in terms of the 
wall-normal velocity v y and vorticity w y and operating on each equation by A , 
results in the following equations [T8j : 



(A.7) 





= A 


Vy 


Wy_ 




Wy_ 



where 
A = 



An 
A21 Ain_ 
-A- 1 Ud x A + A~ 1 U"d x + A- 1 A 2 /R 
-U'd z 





-Ud x + A/R 



(A. 



The term An is the well-known Orr-Sommerfeld operator acting on the wall-normal 
velocity, while A12 and A22 are referred to as the coupling and Squire operators, 
respectively. The eigenvalues of the OS operator can be studied as a function of the 
wavenumber in the a;-direction, which we denote a. 

Appendix B. Spectral methods. 

Differentiation matrices approximate differential operators acting on functions as 
matrices acting on vectors. Finite difference methods track function values at points 
in physical space, while spectral methods approximate functions as finite series of basis 
functions, and keep track of the series coefficients. Because there is some confusion 
in terminology for spectral methods, we briefly review definitions used in this paper. 

In general, spectral methods approximate functions by a truncated series of N 
basis polynomials, Jn(v) = 'Yli 1 j=o a j4 > j- We will refer to formulations with operators 
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that act directly on these polynomials as basis function schemes. Alternatively, spec- 
tral collocation utilizes the fact that there is a one-to-one correspondence between the 
coefficients of that series, a,j, and the values of the function at specially chosen, non- 
uniform grid points, ftf(yj). Each function f(y) can be approximated by its values 
at these special set of points, /at(j/j), and operators are approximated as matrices 
that act on these points. The grid points are chosen so that the error associated with 
the approximation is evanescent, or O {(1/N) N ), which is superior to finite difference 
methods. A particularly clear introduction to spectral collocation is given by Boyd [3]. 

There are two main methods for applying BCs in spectral methods. The first 
is basis recombination, where fields are expanded in a series of basis functions that 
independently satisfy the appropriate BCs [H[6], which are called Galerkin schemes. 
BCs can also be enforced using boundary bordering, which has two variations. The 
first variation enforces homogeneous Dirichlet BCs in spectral collocation schemes. In 
this case, the first and last entries in the vector correspond to the boundary points 
in physical space, and one simply removes the first and last rows and columns of the 
differentiation matrix. Note that this reduces the size of the square differentiation 
matrix by two. A second variation involves using m rows of the matrix to enforce 
m BCs explicitly, which is often referred to as the Tau method. Note that both Tau 
and the Galerkin methods can be implemented with collocation and basis function 
schemes. 

One numerical scheme referred to often in this paper is spectral collocation with 
Chebyshev polynomials. The ij entry of the first derivative matrix, Cheb^j , is given 
I'.v ®- 

( (l + 27V 2 )/6, i=j = 0; 

y I -Xj/ [2(1 - xj)\ , i = j; < j < N; v ' 

{ (-ly+tpi/faixi-xj)], i+j, 

where 

p =PAr=2, Pj = l, je(l...N-l). (B.2) 



The second-derivative spectral collocation approximation is given by the square of the 
first-derivative matrix: Cheb^ » = (Cheb«) . 
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